################################################################
# Replication File: "Backscratching in Banks: Political Cycles in Bank Manager Appointments."
# Author: Jonas Markgraf
#
# Main Analysis contains (in chronological order as appearing in paper):
## -- replication code for all empirical analyses and plots in the paper
#
# R Version: 3.6.0 (2019-04-26)
# platform:  x86_64-apple-darwin15.6.0
################################################################

rm(list = ls())

# load packages
library(eha)
library(survival)
library(flexsurv)
library(texreg)  
library(dplyr)
library(simPH)

# load data

## data for Spanish savings banks
cajas.boards <- read.table(file = "cajas_dta.txt", header = T)  

###############################################################################
############## MAIN RESULTS ###################################################
###############################################################################

# **Table 1: Main Results**:  -----------------------------------------------

## Column 1
fit.main1 <- coxph(Surv(bank.yrsInOffice, bank.event) ~  region.yrsInGov,
  method = "efron",
  data = cajas.boards, subset = bank.chair == 1)
summary(fit.main1)

## Column 2
fit.main2 <- coxph(Surv(bank.yrsInOffice, bank.event) ~  region.postElec,
  method = "efron",
  data = cajas.boards, subset = bank.chair == 1)
summary(fit.main2)

##Column 3
fit.main3 <- coxph(Surv(bank.yrsInOffice, bank.event) ~ region.yrsInGov +
    region.postElec +
    region.coal + 
    bank.polVote +
    bank.roaLag1 +
    as.factor(region.govParty) + as.factor(region.govParty):bank.yrsInOffice +
    frailty(bank.id, distribution = "gaussian"),
  method = "efron",
  data = cajas.boards, subset = bank.chair == 1)
summary(fit.main3)

## Column 4
fit.main4a <- coxph(Surv(bank.yrsInOffice, bank.event) ~ region.postElec +
    region.coal:region.postElec +
    region.coal + 
    bank.polVote +
    bank.roaLag1 +
    as.factor(region.govParty) + as.factor(region.govParty):bank.yrsInOffice +
    frailty(bank.id, distribution = "gaussian"),
  method = "efron",
  data = cajas.boards, subset = bank.chair == 1 & region.govChange == 1)
summary(fit.main4a)

## Column 5
fit.main4b <- coxph(Surv(bank.yrsInOffice, bank.event) ~ region.postElec +
    region.coal:region.postElec +
    region.coal + 
    bank.polVote +
    bank.roaLag1 +
    as.factor(region.govParty) + as.factor(region.govParty):bank.yrsInOffice +
    frailty(bank.id, distribution = "gaussian"),
  method = "efron",
  data = cajas.boards, subset = bank.chair == 1 & region.govChange == 0)
summary(fit.main4b)

## print all results
screenreg(list(fit.main1, fit.main2, fit.main3, fit.main4a, fit.main4b),
  caption = "Cox PH Model", caption.above = T,
  label="T:findingsChair",
  custom.coef.names = c("Years in Government"
    , "Post Election"
    , "Coalition"
    , "Return on Assets$_{t-1}$"
    , "Public Sector Vote Share"
    , "Party: Other", "Party: Socialist"
    , "Post Election * Coalition"),
  custom.model.names = c("Raw Model 1", "Raw Model 2", "Full Model", "After Gov't Change", "No Gov't Change"),
  include.rsquared=F, include.maxrs=F, include.missings=F, include.zph=F,
  reorder.coef = c(1:2,8,3:7),
  stars = c(.01,.05,.1),
  omit.coef = c("bank.yrsInOffice"))

## **Figure 1**: Plot effect of time in government ----------------------------

## simulate values from Table 1, model 3, using 'simPH' package
plot.region.yrsInGov <- coxsimLinear(fit.main3
  , b = "region.yrsInGov"
  , Xj = seq(1, 28, by = 1)) # for year 1--28

## plot
fig1 <- simGG(plot.region.yrsInGov
  , xlab = "Years in Government"
  , lcolour = "black"
  , pcolour = "grey"
  , alpha = .5
  , type = "points")
pdf("figures/effect_govChangeYr.pdf")
fig1
dev.off()

# **Figure 2**: Plot interactions from Table 1 -------------------------------

## prepare coefficients from Table 1, model 4, for plotting

### extract values
fit.main4a$coefficients
relevant.coefs1 <- as.data.frame(coef(fit.main4a)[c(1:2,7)])

### obtain full interaction
relevant.coefs1[3,] <- sum(coef(fit.main4a)[c(1:2, 7)])  

### calculate standard error for interaction
se_interact1 <- sqrt(fit.main4a$var2[1,1] + fit.main4a$var2[2,2] + 2*fit.main4a$var2[1,2])

### extract values for 90% confidence interval
relevant.cis1 <- as.data.frame(confint(fit.main4a, level = .9)[c(1:2,7), c(1:2)]) 

### estimate 90% confidence interval for full interaction
relevant.cis1[3,1] <- relevant.coefs1[3,] - 1.645*se_interact1
relevant.cis1[3,2] <- relevant.coefs1[3,] + 1.645*se_interact1

### bind everything together and clean
exp.int1 <- cbind(relevant.coefs1, relevant.cis1)
colnames(exp.int1)[1] <- "est"

### drop main effect for coalition
exp.int1 <- exp.int1[-2,]
exp.int1

## prepare coefficients from Table 1, model 5, for plotting

fit.main4b$coefficients
relevant.coefs2 <- as.data.frame(coef(fit.main4b)[c(1:2,7)])
relevant.coefs2[3,] <- sum(coef(fit.main4b)[c(1:2, 7)])
se_interact2 <- sqrt(fit.main4b$var2[1,1] + fit.main4b$var2[2,2] + 2*fit.main4b$var2[1,2])
relevant.cis2 <- as.data.frame(confint(fit.main4b, level = .9)[c(1:2,7), c(1:2)])
relevant.cis2[3,1] <- relevant.coefs2[3,] - 1.645*se_interact2
relevant.cis2[3,2] <- relevant.coefs2[3,] + 1.645*se_interact2
exp.int2 <- cbind(relevant.coefs2, relevant.cis2)
colnames(exp.int2)[1] <- "est"
exp.int2 <- exp.int2[-2,]
exp.int2

## plot the obtained estimates for Table 1, model 4 and 5

pdf("figures/interaction_log-hazardratio.pdf")
par (mar=c(4,6,0,0))
plot (c(1,3), c(floor(min(exp.int1[2])), ceiling(max(exp.int1[3]))), type="n", bty="n"
  , ylab=""
  , xlab=""
  , axes=F)
y_axis <- seq(floor(min(exp.int1[2])), ceiling(max(exp.int1[3])), by=1)
axis (2, at = y_axis, labels = parse(text = paste("e^", y_axis)), las = 1)
mtext (text="Relative Hazard", side=2, line=2.5, cex = 2)
mtext (text=c("Single Party", "Coalition"), side=1, at=c(1.5,2.5), line=0.5, cex = 2)
abline(h=0, lty =2)
points (xy.coords(c(1.4,2.4), exp.int1[,1])
  , pch=19, col=c("black","black"), cex=2)
segments (x0=c(1.4,2.4), x1=c(1.4,2.4), y0=exp.int1[,2], y1=exp.int1[,3]
  , lwd=3, col=c("black", "black"))
points (xy.coords(c(1.6,2.6), exp.int2[,1])
  , pch=19, col=c("grey","grey"), cex=2)
segments (x0=c(1.6,2.6), x1=c(1.6,2.6), y0=exp.int2[,2], y1=exp.int2[,3]
  , lwd=3, col=c("grey", "grey"))
legend("topright", legend = c("New Government", "Reelected Government"), col = c("black", "grey"), lty = 1
  , pch = c(19, 19), box.lty=0, cex = 1.5)
dev.off()
